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ABSTRACT 


COUPLED BENDING-TORSION STEADY-STATE RESPONSE OF PRETWISTED, 
NONUNIFORM ROTATING BEAMS USING A TRANSFER-MATRIX METHOD 

Using the Newtonian method, the equations of motion are developed 
for the coupled bending-torsion steady-state response of beams 
rotating at constant angular velocity in a fixed plane. The resulting 
equations are valid to first order strain-displacement relationships 
for a long beam with all other nonlinear terms retained. In addition, 
the equations are valid for beams with the mass centroidal axis offset 
(eccentric) from the elastic axis, nonuniform mass and section proper- 
ties, and variable twist. The solution of these coupled, nonlinear, 
nonhomogeneous, differential equations is obtained by modifying a 
Hunter linear second-order transfer-matrix solution procedure to solve 
the nonlinear differential equations and programing the solution for a 
desk top personal computer. The modified transfer-matrix method was 
verified by comparing the solution for a rotating beam with a 
geometric, nonlinear, finite-element computer code solution; and for a 
simple rotating beam problem, the modified method demonstrated a 
significant advantage over the finite-element solution in accuracy, 
ease of solution, and actual computer processing time required to 
effect a solution. 
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CHAPTER 1 


INTRODUCTION 

At the Langley Research Center there are several wind tunnels 
which are driven by rotating propeller or fan blade systems. Nearly 
all of these tunnels have been in operation for several years; and in 
order to ensure safe and continued operation, the existing blades will 
require a thorough evaluation of the steady-state response under the 
current operating conditions. To determine the steady-state response 
of rotating blades, effective and efficient analysis techniques are 
required that consider the unique features common to rotating systems. 

An important consideration in the structural analysis of a 
rotating beam is the effect of the coupling between the centrifugal 
force and the elastic stiffness (deflections) of the beam. This 
effect provides a stiffer structural member in both the beamwise and 
chordwise directions. To solve this problem and take advantage of the 
centrifugal -stiffening effect is difficult due to the necessity of 
having to know the deformed state of the beam in order to apply the 
steady-state centrifugal loads. 

The literature is extensive with current and classical papers 
concentrating on the rotating beam problem. Among the more notable is 
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the linear analyses of Houbolt and Brooks [1]*. In their paper, they 
systematically develop the governing differential equations for the 
coupled bending and torsion of pretwisted blades assuming shear and 
rotary inertia are negligible. The warping rigidity for a slender 
beam was implicitly considered by incorporating a St. Venant torsion 
term in the internal elastic torque equation. Hunter [2] developed a 
system of linear coupled differential equations for small displace- 
ments where the axial displacement and centrifugal force equations 
easily uncouple, and allows the axial component equations to be solved 
separately from the transverse component equations. These equations 
illustrate the basic but essential physics of the rotating beam and 
were developed specifically to demonstrate an integrating-matrix 
method for solving homogeneous differential equations. Shear and 
rotatory inertia effect were considered by Stafford and Giurgiutiu [3] 
for the uncoupled uniform vibrating beam rotating in a fixed plane. 
They found that these effects were more noticeable for blades 
operating at higher speeds (turbomachinery blades). These and other 
similar papers, however, center around predicting the natural vibra- 
tion characteristics [2,3,4] of rotor blades (rotating beams). The 
concern with the natural vibration problem is to accurately determine 
mode shapes and natural frequencies of rotors to avoid resonance. The 
natural vibration problem, in effect, solves the eigenvalue problem 
resulting from eliminating the nonhomogeneous part of the governing 
equations of motion. Solutions of the eigenvalue problem are 
important in the design of rotating beams, but solutions for the 


* The numbers in brackets indicate reference. 
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steady-state response of the rotating beam are also important for 
determining strength for fatigue evaluation and distortions for 
clearance evaluations. 

The purpose of this thesis is to address the development and 
solution of the governing equations of pretwisted rotating beams under 
the influence of steady-state aerodynamic and inertial loadings (the 
nonlinear nonhomogeneous boundary-value problem) using a transfer- 
matrix method. The transfer-matrix solution method has been selected 
because of its capability to account for the variation of the beam's 
geometry (stiffness and mass) and loading along the span of the beam. 

Chapter 2 first develops the kinematics of the longitudinal 
deformations of the beam including uniform extension, bending about 
two axes, and uniform twist. Second, the stress resultants for a 
linear elastic material are developed using the longitudinal strain. 
Following the development of the internal stress resultants, the 
acceleration of a differential element where the elastic and 
centroidal axes are not coincident is derived for a steady-state 
rotating element. By application of D'Alembert's principle to the 
resulting steady-state inertial loads, the equilibrium equations for a 
differential element under transverse and torque loadings are 
developed. Chapter 3 formulates the resulting 12, coupled, nonlinear, 
first-order, ordinary differential equations derived in Chapter 2 into 
a classical boundary-value problem. Since the intended application of 
this development is to provide a complete analytical tool, Chapter 4 
presents a nonlinear solution method for solving the derived 
equations, and Chapter 5 discusses the computer programs necessary to 
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implement the solution algorithm. Chapter 6 describes an application 
of the theory, solution algorithm, and computer programs necessary to 
analyze a new fan blade design for the 7- by 10-Foot Wind Tunnel at 
the Langley Research Center. Results from the transfer-matrix approach 
are evaluated by comparison with finite-element results. In Appendix 
A, a derivation of the second-order Hunter transfer-matrix method is 
presented, and Appendix B gives the details of calculating the 
necessary cross-sectional integrals for a typical 7- by 10-Foot Wind 
Tunnel fan blade. 


4 



CHAPTER 2 


MATHEMATICAL FORMULATION 

The system of displacements for a generalized flexure theory is 
derived based upon the assumption that regardless of the loading, the 
original shape of the cross section is unaltered during deformations. 
Thus, the geometric dimensions of every plane normal to the longitudi- 
nal axis of the beam remains unchanged. Note, however, that this 
assumption precludes any type of out-of-plane cross-sectional defor- 
mations (neglect warping and shear deformations). Following from these 
displacements, the longitudinal strains are developed as is done by 
Houbolt and Brooks [1] for a beam under both lateral, extensional, and 
twisting motion. After expressing the inertial loads on a beam in 
terms of the deformations, the derivation then considers the 
equilibrium of a deformed beam in terms of its geometry and resultant 
loads, thus, yielding a set of nonlinear governing differential 
equations for the longitudinal, transverse, and torsional deformations 
of an elastic structure. 


2.1 Analysis of Displacements 

The cross-sectional plane defined by x equals a constant (see 
figure 2.1) is allowed to undergo u, v, and w displacements of the 
elastic axis. In addition, the cross section is allowed to rotate ©y 
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and 0 Z about the y and z coordinate axes, respectively. The cross- 
sectional twist, <f> , is measured about the elastic axis, x. By passing 
a plane through the beam that is normal to the cross-section's x-axis, 
the location of a point, P, on the pretwisted cross-section may be 
described for both the undeformed (figure 2.2(a)) and deformed, 

(figure 2.2(b)) beam. Attached to the cross-section is a (n,s) 
coordinate system that moves with the deforming cross-section and is 
located at the shear center. The angle 0 locates the positive n-axis 
which locates the principal minor axis of inertia relative to the 
positive y-axis. 

The position vector r to the point P in the y-z and n-c 
coordinate systems is given by 



r > 
X 


1 

0 


0 

— 


r > 
X 

r = 

n 

> = 

0 

COS 

(B) 

sin 

(0) 

i 
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c 

V* J 


0 

-sin 

(0) 

cos 
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Z 

^ j 
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P 
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0 
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r 

X 

r = < 

y 

► = 

0 

cos 

(0) 

-sin 

(0) 
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n > 
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0 

sin 

(0) 

cos 

(0L 


c 


v. J 


The rate of change, using primes to denote differentiation with 
respect to x, of the point P with respect to the longitudinal axis is 
then given by 


y 


► 

l 

J 



r~sin ($) 

-cos (ar 


fnl 

0' 

COS (0) 

-sin (3L 

| < 

[J 


(2.3) 
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Using equation (2.1) in (2.3) yields 



► = S' < 


b'J 


bJ 


Next by imposing the u, v, and w displacements at the elastic axis and 
also rotating the cross section by the angle <j>, the position vector 
£1 (see figure 2.2) is expressed in terms of the original position 
vector r and the cross-sectional displacements. Recall that the 
fundamental assumption of this development is that the original cross- 
section does not change its shaper this implies that the vector £ 
relative to the n-c system is unaffected. Thus, only the kinematics 
of the cross section in the y-z system are evaluated. 


The xj component of the £j vector is composed of contributions 
due to extension and bending about the y and z axes. Therefore, by 
superimposing the extension (figure 2.3(a)) and bending (figures 
2.3(b) and 2.3(c)) effects, the xi component is 


X 1 


x + u 




(2.5) 


or by using primes to denote differentiation with respect to x, 

xi = x + u - y v 1 - z w' (2.5) 

The other two components of £j, y^ and zj_, are found by using 
equation (2.2) with the pretwist angle, 0, replaced by the total 
twist, (0+<|>). 
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cos(p + <|>) -sin(p + <|>) 
_sin(p + 4 >) cos(p + 4>) _ 


( 2 . 6 ) 




In addition, if the usual small -angle assumptions are made for 4> and 
using double-angle trigonometric identities, then equation (2.6) 
becomes 


M 


[cos (p) - <|> sin (p)] [- sin (p) + 4> cos ( e ) ] 

UJ 

► = 

_[sin (p) + 4> cos (p)] [cos (p) - <j> sin (p)] _ 



(2.7) 


Expanding equation (2.7) and using equation (2.1) leads to 



( 2 . 8 ) 


Superimposing the v and w displacements of the elastic axis and 
including equation (2.5), the position vector r± becomes 


r y 

XI 


^ i i ^ 

x + u -y v -z w 

<yi 


v + y -z <|> 

z i 

V. J 


w + z +y 4) 


(2.9) 


Using equation (2.9), the rate of change of £i with respect to x is 
given by 







H-* 

+ 

c 

1 

<< 

< 

1 

<< 

< 

1 

Nl 
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1 

N 

s 

11 ' = « 

yi' 

ii^ 

v'+y'-z' 4>-z4)' > 


z i’ 

C J 


w' + z 1 + y 1 4> + y 4>' 

^ J 


( 2 . 10 ) 
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Collecting terms and using equation (2.4), equation (2.10) becomes 


II 


* — 


0 

xi 


r l + u 1 - y( v"+ S' w 1 ) - z( w 1 1 - e' v ' ) 

n' 

> = < 

v' - y g' <|> -z(e'+ ' ) 

zi' 

j 


w' + y ( 6 ' + 4» ' ) - z e' <|> 

L J 


) ( 2 . 11 ) 


Thus, after deformation, the position vector of a point, P, is given 
by equation (2.9), and its derivative with respect to x is given by 
equation (2.11). 


2.2 Analysis of Strain 

From the displacement analysis results, the longitudinal strain, 
e, that is developed in a longitudinal fiber, is derived analytically. 
Let the point P under displacements (u,v,w) of point o (elastic axis, 
see figure 2.2(b)) and rotations ( 4>, v', w') of the cross-section, be 
displaced to P' (xi, yj, z^) on the rotated cross-sectional plane as 
given by equation (2.9). Then the differential length, dsi, 
correponds to 

(dsx) 2 = (dxi ) 2 + ( dy i ) 2 + (dzi) 2 (2.12) 


or 

(dsi) 2 = (xi ') 2 + (yi’) 2 + (zi 1 ) 2 (2.13) 

Performing the indicated algebra and using equation (2.11), each of 
the components on the right hand side of equation (2.13) has the 
following form 


9 



(X! 1 ) 2 = 1 + 2u' + (u ' ) 2 
+ (- 2v ' 2w' 3 2u ' v ‘ ' - 23 'u'w 1 )y 
+ (- 2w''+ 2v ' e ' - 2u'w l, + 2 3'u'v')z 
+ (2v"w" - 23'vV - 2( 3 1 ) 2 v'w' + 2b'w ' ) yz 
+ ( (v") 2 + 23Vv" + (3') 2 (w') 2 )y 2 

+ ( (w") 2 - 23'v'w" + (3') 2 (v') 2 )z 2 (2.14) 

(yi ' ) 2 = (v') 2 + (-2 3 V<t>)y + (-2 3'v' - 2 <t>' v' ) z 
+ (2(3 1 ) 2 4> + 2 (3'<p'<t>)yz + ( (3 1 ) 2 (4>) 2 ) y 2 
+ ( ( 3 1 ) 2 + 2 3V + (V) 2 )z 2 (2.15) 

(zi')2 = (w 1 ) 2 + ( 23'W* + 2w ' <(> ' )y + ( - 2 B'wVz 

+ ( -2(3' ) 2 <j> - 2 3V<I> )yz + ( (3') 2 + 23* 4>‘ + V) 2 )y 2 
+ ( (3 1 ) 2 ( ) 2 )z 2 (2.16) 

Substituting equations (2.14)-(2.16) into equation (2.13) leads to 

(si') 2 = 1_ + 2u' + (u ' ) 2 + (v 1 ) 2 + (w 1 ) 2 
+ ( -2v ' 1 - 2u'v'' - 2u , w'3 l - 2 3'v'<f> + 2<j>V )y 
+ ( -2w' 1 - 2u ’ w ' ’ + 2uV3' - 2 3 1 w ' 4 > - 2<f>V )z 
+ (2v ,, w" - 2v' 3 V 1 + 2w , 3 , w" - 2v‘w'(3') 2 )yz 
+ ((v") 2 + 2w 1 3 1 v ' ' + (w ' ) 2 ( 3 ' ) 2 + ( 3' ) 2 + 23V 
+ V) 2 + (3 1 ) 2 (<t>) 2 )y 2 

+ ((w") 2 - 2v'3 , w ,, + (v') 2 (3') 2 + ( 3 1 ) 2 + 23 V 
+ V) 2 + ( 3 1 ) 2 (<j> ) 2 )z 2 

(2.17) 
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For this development, the assumption is that the problem formulation 
be restricted to small displacements. Thus, the second-order terms 
in equation (2.17) are much less than unity; therefore, only the 
first-order displacement terms, (underlined in equation 2.17), are 
retained. 

^ = [1 + 2u ' -2y v -2zw" + (y 2 + z 2 )((g') 2 + 2e'<J>')] 1/2 

(2.18) 

Taking the initial length of ds to correspond tou = v= w = (j> = 0, 
equation (2.18) reduces to 

^ = [1 + (y 2 + z 2 ) (e*) 2 ] 1/2 (2.19) 

dx J 

Making use of the fundamental definition of longitudinal strain to the 
first order in dx 


dsi - ds 
3s 


( 2 . 20 ) 


and equations (2.18) and (2.19), the following result is obtained 

e = { 1 +2u ‘ -2yv“ -2zw' '+(y 2 + z 2 )(g' 2 +26 ' <j>) ) • 5 - {1 +(y 2 +z 2 )e' 2 } - 5 

(1 + (y 2 + z 2 ) (3 1 ) 2 ) - 5 

( 2 . 21 ) 

Noting that each of the square root terms in equation (2.11) is of the 
form 


(1 + 6 ) 1/2 


( 2 . 22 ) 
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which can be, for small 6, approximated by 


1 + 6/2 (2.23) 

Making use of equation (2.23), then equation (2.21) can be reduced to 

u 1 - y v 1 1 -zw 1 1 + (y 2 + z 2 ) 6 1 <j> 1 

1 + (y 2 + z 2 ) i|^ 2 (2.24) 

A closer evaluation of the second term in the denominator of equation 
(2.24), for most rotating beams that are used as fan blades and 
propellers, reveals that the (e 1 ) 2 term is usually on the order of 
1° per inch. This can, for most geometries, reduce the second term 
to the order of 0.01 which is much less than unity; therefore, 
it is neglected. This yields the final linear expression for the 
longitudinal strain of 

e = u' -yv” -zw” + (y 2 + z 2 )e'<|>' (2.25) 

As pointed out in reference [1], the above expression for the 
longitudinal strain includes the effects of four general types of 
motion: uniform longitudinal extension of the cross section; rotation 

of the cross section due to two axes bending; and rotation of cross- 
sectional planes relative to each other about the elastic axis due to 
twisting of the beam. 

2.3 Analysis of Inertial Loads 

The acceleration of a rotating beam is determined at the centroid 
of the segment, see figure 2.4. For this derivation, the axis of 
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rotation is assumed to coincide with the z-axis, which is normal to 
the elastic axis. Thus, the location of the center of mass on the 
cross section is given by 

Jjn = Xmi* W + ^n!i < 2 - 26 > 

where the components of r„,, in the deformed configuration, are 
determined by using equation (2.9). 

^m = x m + u " ^m v ' “ z m w ' (2.27) 

\ “ ^m + v “ z m ♦ (2.28) 

?m = z m + w + ym * (2.29) 

The origin of the elastic axis is assumed to coincide with an 
inertial frame such that the elastic-axis coordinate system is seen as 
just a rotating system with respect to this inertial frame. Thus, the 
position vector, equation (2.26), describes the location of the center 
of mass of each cross-section explicitly in terms of the undeformed 
location and the displacement components; therefore, the time deriva- 
tives (denoted by dots) of the position vector are easily determined. 
The velocity of the center of mass is given by 

• • • • • • • 

rm = *m i + X m + Y m j + Y m j + Z m k + Z m k (2.30) 


13 



Since the U j, and k system is just rotating with respect to the 
inertial frame, the time derivatives of the bases vectors are 


i_ = n j 

£ = -fl 1 

k = 0. 


(2.31) 

(2.32) 

(2.33) 


Also, the time derivatives of the components of the position vector rjp 
are 


^m = u - ym V - z m w' 

Y m " v - z m i 

• • • 

Zm = w + ^m 

Using equations (2.31 )— (2.33) , equation (2.30) becomes 


(2.34) 

(2.35) 

(2.36) 


Im “ f^m “ "^m ) 1 + ( ^m + n ^m ) 1 + ^m Ji (2.37) 

Following the same approach, the acceleration of the center of mass of 
a cross section, assuming a constant rotation, is 


Im ~ Z Q T m 


- + (T m + 2 si \ - Sl 2 Y m )j + Z m k (2.38) 
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where 


x m u 




v - 


z m w 


(2.39) 


Y 


v - z m ♦ 


(2.40) 


Z n = w + y m * (2.41) 

Upon substituting equations ( 2 . 39 ) - ( 2 . 41 ) into equation (2.38), the 
acceleration of the center of mass becomes 


£m = Cu " ymv' - W* - 2 o(v- z m J) - fl 2 (x m +u-y m v' - z m w')] i 

[v - z m $ + 2 n(u -y m v' - Z m w') - n 2 (y m + v - z m <|>) ] j 

[ w + y m i ] l< (2.42) 

For the purpose of this analysis, the motion of concern is the steady- 
state response of the beam. Thus, all time derivatives are zero; and 
equation (2.42) becomes 

£m = -a 2 (x m + u - y m v'- z m w') - n 2 (y m + v - z m *) j 

= A x i + Ay j (2.43) 
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2.4 Development of Governing Equation 


2.4.1 Derivation of Equilibrium Condition 
For the derivation of the equilibrium conditions of a rotating 
beam as shown in figure 2.5, attention is restricted to the steady - 
state response of a beam that is rotating at a constant angular 
velocity under steady-state aerodynamic lift (P z ), drag (Py), and 
torque (q x ). The aerodynamic loads are assumed to be resolved to the 
elastic axis as varying distributive loads. 

A differential beam segment is shown in figure 2.6 with its 
associated loads: internal, inertial, and aerodynamic. For the 

steady response of a rotating beam, the differential inertial loads 
are assumed to act at the center of mass of the differential segment. 
Using D'Alembert's principle of dynamic equilibrium, the inertial 
loads are treated as if they were applied static forces, but in the 
opposite sense. Referring to figure 2.6(a), the axial equilibrium 
condition is 


T+ — dx-T - P A x dx = 0 (2.44) 

dx 

Upon substituting the ji-component of equation (2.43) into equation 
(2.44) and simplifying, the governing equation for the tension becomes 

— + p fi 2 (x m + u - y m v' - z m w') =0 (2.45) 

dx 

Summing the forces in the y-di recti on, figure 2.6(b), the following 
equilibrium equation is derived 
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Using equation (2.43), equation (2.46) becomes 


riV 

^ + p Q 2 (v + y m -z m 4 ) + Py = 0 (2.47) 

Similarly, the z-force equilibrium equation is derived, see figure 
2.6(c). 


dV, 

V 2 — dx - ^2 ^ = ^ 


(2.48) 


or 




(2.49) 


The sum of the moments in the x-y plane, about the center of the 
differential element, figure 2.6(b), yields the following 


dM 


dx 


dVv 


dx 


M 2 + dx - M z - pA x (y m -z m <j>)dx + — (Vy + dx) + — Vy 


-v 1 ^ (T + dx) - v'£ T - 0 
2 dx 2 


(2.50) 


Dividing through by dx, using equation (2.43) for A x , and taking the 
limit as dx goes to zero, equation (2.50) reduces to 



(2.51) 


dM 7 o * o i 

-p nZ ( x m ^m + u ^m “ v ^m 2 " w 2 m y m " *m z m 4> ~ u z m^ 

+ v 'y m 2 m ♦ + w z m 2< f > ) + V y - v’T = 0 
Similarly the moment equilibrium equation in the x-z plane is derived 

(nr + p a2(x i» z m + u z m -*'*■ z nT w ' z m 2 *>.%**“ y m * 

-v'ym 2 * -w'ym - V 2 t w'T = 0 (2.52) 


From figure 2.6(a), the equilibrium equation for the y-z plane 
(torque) is 

~ + dx “ p ( z m v + ^m z m “ z m 2 4> + ym v 4> + ^m 2 4 > - z m y m 4>^) = 0 

(2.53) 

There are now 6 equations, (2.45,47,49,51-53), and 10 unknowns (T, V y , 
V z » M x » M y , M z , u, v, w, <|)) which indicates that additional relation- 
ships are required. The additional equations can be developed from 
the internal elastic equilibrium (stress resultant) constraint imposed 
on a cross section that the integrated stresses balance the resultant 
loads at the elastic axis. 

2.4.2 Internal Elastic Loads 

By assuming the longitudinal material orientation of the rotating 
beam behaves as a linear elastic material, the constituent relation- 
ship can be described by Hooke's law. 

o = E e (2.54) 
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where E is the apparent longitudinal modulus of elasticity, and e is 
given in equation (2.25). 

From the above definition, the stress distribution over the 
cross-section may be resolved into the internal resisting loads at the 
elastic axis, see figure 2.7. 

The axial, T, and bending. My and M z , internal loads are given by 



(2.55) 



z dA 


(2.56) 


M z «- f o y dA 


(2.57) 


By introducing the minus sign on the M z component, a tensile 
longitudinal stress will be produced when a negative M z is introduced; 
this makes the sign of the normal stress consistent with the sense of 
the internal load, see figure 2.6a. As pointed out by Houbolt and 
Brooks[l], the selection of the elastic axis as the primary reference 
axis was to eliminate the shearing stresses that are associated with 
the longitudinal stresses from contributing to the total resisting 
torque, M x . Because of the pretwist, the longitudinal stress has a 
component that contributes to the internal torque (see figure (2.8)). 
From figure (2.8), it is seen that the inplane components of the 
normal stress produce a torque 
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M xn = f a + ^) ( 3 ' + 4 .') dA 


(2.58) 


Substituting equation (2.1) into equation (2.58) and noting that the 
Jocobian relating the differential areas in the y-z and n-c systems is 
unity, equation (2.58) becomes 

M xn = fa (S' + ♦')[ y 2 + z 2 ]dA (2,59) 

As is done in [1], equation (2.59) is combined with the St. Venant 
twisting which leads to the following equation for the total resisting 
internal torque. 


M x = fa (e 1 + <f>')[ y 2 + z 2 ]dA + GJ 4 > ' 


(2.60) 


For convenience, the internal resisting loads are restated again using 
equation (2.25) 


T = E f ( u ' -yv" -zw' 1 +(y 2 + z 2 )e'<f>')dA ' 2 ‘ 61 ' 

My = E f z(u' -yv" -zw" +(y 2 + z 2 ) 0 1 ' )dA ^ 2 * 62 ^ 

M z =-E f y ( u ' -yv" -zw" +(y 2 +z 2 ) e’ <j>' )dA (2.63) 
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M x - E f (y 2 + z 2 )(g - 

+ 4> ' ) ( u ' -yv" - 2 W 1 ' +(y 2 + z 2 )8'<f>' 

) dA + GJ<f> ' 
(2.64) 

By making use of the following 

definitions, 


> 

II 

<-> 

dA 

(typical units in 2 ) 

(2.65) 

% 

II 

*-> 

z 2 dA 

(typical units in 4 ) 

(2.66) 

l zz = / 

y 2 dA 

(typical units in 4 ) 

(2.67) 

H— 1 
<< 
Nl 

II 

yz dA 

(typical units in 4 ) 

(2.68) 

’p = S 

(y 2 + 

z 2 ) dA (typical units in 4 ) 

(2.69) 

p i ■ / 

y dA 

(typical units in 3 ) 

(2.70) 

S 

z dA 

(typical units in 3 ) 

(2.71) 

p 3 = r 

j 

z(y 2 + 

z 2 ) dA (typical units in 5 ) 

(2.72) 
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P4 = J* y(y 2 + z 2 ) dA (typical units in 5 ) 


(2.73) 


P 5 = J 4 (y 2 + z 2 ) 2 dA (typical units in 6 ) (2.74) 


the integrations, as indicated, in equations (2.61)-(2.64) can be 
performed to yield the following 


T = E [ A u' - Pi v" 

- P 2 W 1 ' 

+ ip 3 1 <t>' ] 

(2.75) 

My = E [ P 2 u‘- Iy Z v’ 

1 - I w 1 

A yy w 

' + P 3 6 ' ♦' ] 

(2.76) 

M z = E [-Pi u'+ I 2Z v 1 

' + I w' 

Xyz W 

" - p 4 8 1 ♦' ] 

(2.77) 


M x = GJ <>' + E [ I p 3 1 u * - P 4 8 ' v" - P 3 e' w" + P 5 ( 0')2 
+ Ip 4 > 1 u 1 - P 4 <}> 1 v ' 1 - P 3 <t> 1 w' 1 + P 5 ( 4> ' ) ^ ] 

(2.78) 

Equations (2.75)-(2.78) define the internal elastic equilibrium which 
relates the internal forces to displacements. This now increases the 
number of equations to 10 , and thus, ensures a compatible system of 
equations and unknowns. 
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CHAPTER 3 


FORMULATION OF GOVERNING EQUATIONS AS A CLASSICAL 
BOUNDARY-VALUE PROBLEM 

The preceding analysis has furnished a set of coupled nonlinear 
differential equations. These equations together with constraints or 
conditions on the boundaries comprise the classical boundary-value 
problem. This chapter is an organized summary of the previously 
derived governing equations recast and formatted in a form that yields 
a system of first-order differential equations. 


3.1 Differential Equations 

To summarize, the six equilibrium equations (2.45), (2.47), 
(2.48), (2.51), (2.52), and (2.53) are restated. 


dT = 
dx 

dVy _ 

IT ~ 


dV z . 
W " 


p n 2 (x m + u - y m v 1 - 
- P fl 2 (v + y m - z m <f>) 

- Pz 


-m 


w') 



(3.1) 

(3.2) 

(3.3) 


dMx 

dx 


' q x + p (zm v + y m z m 


z m 2 ♦ + y m v ♦ + y m 2 - z m y m *2) 


(3.4) 
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(3.5) 


dM v o . i o 

= * p z m + u z m “ v ^m z m " w z m 2 + x m ^m ♦ + Wm ♦ 

- v'y m 2 4 > - w'ym z m ♦ > + v z ’ W ' T 

a p ^ 2 ( x m ^m + u ^m " v -^m 2 “ wz m " x m z m 4 " u z m ♦ 

+ v 'ym z m ♦ + w ' z m 2 ♦) “ v y + V ' T (3 * 6) 

Similarly, the stress resultant equilibrium equations (2.75-78) are 
restated. 


T 

= E [ A u 1 - Pi v 1 1 

- P 2 w" 

+ I p 6' <t>' ] 

(3.7) 

M y 

■ E [ P 2 u'- I yz v' 

' - ! yy w ' 

" + P 3 6' 4> ' ] 

(3.8) 

M z 

- E [-Pi u 1 + I zz v 1 

' + Iyz W 

l — 1 

ca 

CL 

1 

(3.9) 


M x = GJ </>' + E [I p 3' u' - P 4 B 1 v" - P 3 S' w" + P 5 (6 1 ) 2 4 >' 

+ Ip <t> 1 u 1 - P 4 <(> 1 v' 1 - P 3 4 > 1 w 11 + P 5 (<f>') 2 ] 

(3.10) 


Since the proposed solution method utilizes matrix methods, equations 
(3.7-10) are written as first order differential equations in matrix 
form by letting 


M 

1] 



f 

O 


Iw'J 


rH 

1 

1 

0 

1 



(3.11) 
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" EA 

e e' i p 

E P 2 

-E 

P 1 " 


V" 


"T ' 

E P'lp 

GJ + C<|> 

E 3 ' P 3 

-E 

3 ' P4 


♦ ' 


M x 






i 


>= < 

> 

E P 2 

E 3 ' P 3 

E I yy 

-E 

Vz 


V 


| M y 

--E Pi 

-E 3 ' P 4 

-E I yz 

E 

hi - 




S. 


(3.12) 


where Ccj> = E Ip u' - E P4 0 Z ' + E P3 Gy 1 + E P5 4>'+ E P5 ( B ' ) ^ 
Equation (3.12) is of the form 


V " 


— 1 

♦ ' 

< 

II 



M x 

> 

V 


»y 



M Z 

^ Z J 


(3.13) 


Since equation (3.12) is nonsingular, the solution of (3.12) or (3.13) 
can be obtained. 


V > 



"T " 



r 1 

-1 

M x 


11 

o| 

< 

> 

Q y' 


L J 


M y 




f 

2 

N 


(3.14) 
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Figure 3.1 is a summary of the twelve first order differential 
equations in matrix form. These 12 equations are now in a form that 
is easily integrated by a transfer-matrix method. The general 
character of the equations is that four of the equations result in 
solutions for the unknown displacement (u, v, w, <|>), and six of the 
equations yield solutions for the generalized loads while the 
remaining two equations (3.11) are definitions to reduce the second 
order differential equations (3.7-10) to first-order differential 
equations. 


3.2 Boundary Conditions 

As is the case with most boundary-value problems, the types of 
boundary conditions that can be specified are generally either 
geometric, natural, or a complicated combination of both. One special 
feature of the transfer-matrix solution method is that it can easily 
handle any combination of homogeneous or nonhomogeneous, natural or 
geometric boundary conditions. For the rotating beam, the boundary 
conditions are geometric at the fixed support end and natural at the 
free or cantilevered end. 

In terms of the coordinate system shown in figure 2.7, the 
boundary conditions for complete fixity (constrained against all 
displacements) at the support are 

u = v = w = <j> = Oy = 0 Z = 0 (3.15) 

For a free end ( free to displace) the boundary conditions are 

T = V y = V z = M x = M y = M z = 0 (3.16). 
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If there were, for example, a tip mass attached to the rotating beam, 
the axial force as well as the shear forces at the free end may be 
nonzero which results in nonhomogeneous conditions. Again, this 
apparent complication, as is shown in the applications sections, is no 
more difficult to handle than the homogeneous conditions. 

The above boundary conditions can also be represented in a matrix 
format. 


Fixed end (geometric) - 
“1 00000000000 
010000000000 
001000000000 
000100000000 
000010000000 
000001000000 


fu^l 


V 

w 



Mx 

My 


[ 0 ] 


( 3 . 17 ) 


27 



Free end (natural ) - 


0 00000100000 
000000010000 
000000001000 


000000000100 


00000000000 


28 


CHAPTER 4 


TRANSFER-MATRIX SOLUTION METHOD 

A method is developed for numerically computing the solution of a 
set of coupled, first-order, nonhomogeneous, nonlinear differential 
equations. The solution technique employs an integration similar to 
the Hunter method* which is extended here to a fourth-order 
Runge-Kutta scheme to improve the numerical accuracy when limited 
station properties are available. 

The technique of using the transfer matrix to solve a wide class 
of boundary-value problems can readily be found in the literature. 
However, the earlier development of the transfer matrix was primarily 
employed to solve the homogeneous problem [5,6]. This approach proved 
to be a useful analytical tool in determining the natural modes and 
frequencies of vibrating beams, wings, and propellers. Other matrix 
methods such as the integrating-matrix [2] have also been employed, 
with a great deal of success, to solve the nonhomogeneous problem. 

More recent application of the transfer-matrix method [7], addresses 

* The Hunter method, which was developed in 1974, is an 
undocumented transfer-matrix solution process that can easily account 
for the nonhomogeneous part of a set of first-order differential 
equation. A complete development of the method is presented in 
appendix A. 
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both the nonhomogeneous and the homogeneous boundary-value problem. 
These methods, however, require an extensive library of transfer and 
point matrices as well as some manipulation of the matrices to satisfy 
and impose boundary conditions. In addition, each new class of 
boundary-value problems requires a special technique to derive the 
transfer-matrix from the differential equations. In the present paper, 
the presented method allows for a more direct solution of the 
boundary-value problem by formulating the transfer matrix directly 
from the differential equations and solving the linear equations 
completely in two passes over the solution domain. The first pass is 
required to satisfy all boundary conditions, both natural and 
geometric; and the second pass computes the solution at the selected 
stations within the closed interval defined by the boundary points. 

4.1 Fourth-order Runge-Kutta Transfer Matrix 
The Hunter transfer-matrix (appendix A) method can be shown to 
parallel a second-order Runge-Kutta integration. Using this fact, an 
improvement on the accuracy of the Hunter method is easily obtained by 
directly applying a fourth-order Runge-Kutta to equation (A.l). 

{Y'} = [A] {Y} + {B} (4.1) 

Rewriting equation (4.1) in the following functional form 

{ Y 1 } = F (x, Y) (4.2) 

and using a standard fourth-order Runge-Kutta integration [8], the 
value of {Y} at the (i+1) station in terms of the value of {Y} at the 
i station is 
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{ Y-j +i> = {Yj} + 1/6 ({bj> + 2 {b 2 } + 2 (b 3 ) + (b 4 >) (4.3) 

where the { b} 1 s are defined to be 

(bj) = hi F ( x-j , { Y i } ) (4.4) 

{b 2 > = hi F(Xi + 1/2 h i} (Yi } + 1/2 {bj}) (4.5) 

{b 3 } = hi F(Xi + 1/2 hi , (Yi > + 1/2 {b 2 }) (4.6) 

{b 4 } = hi F(xi + hi, {Yi } + (b 3 >) (4.7) 

Using equation (4.2) in equations (4.4-7) yields 

{b!> = ^ [Ai] {Y f } + h f (Bi> (4.8) 

h * ^ 

{ b 2 ) = hi [J] {Yi> + -L- DG [Ail {Yi } 
hi ^ 

+ hi m + ~ OT {Bi> (4.9) 


2 3 

{b 3 } = hi [Ai] (Yi> +-1- [A] [A] {Ti } +-j- [A] [A] [Ai] (Y^ 

2 3 

+ hi CB} + ^- [A] {¥} + -^- [A] [A] (Bi > (4.10) 

h-3 

{b 4 } = hi [Ai+j] { Yi } + hi 2 [Ai+j] [A] { Y n - } + - — [A n * ] [A] [A] { Y n - } 

hi 4 

+ — — [Ai+j] [A] [A] [Ai] (Yi } 

+ hi tBi+l) + "i 2 Ui+O m + ^ [A 1+1 ] [A] <B) 

hi* 4 

+ -4- CAi+i3 [A] m {Bi> 

4 ( 4 . 11 ) 
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where the [X] and {!} matrices are the average of the [A] and {B} 
matrices at the i th and i+1 station, respectively. 


Substituting equations (4.8-11) into equation (4.3), collecting 
terms, and developing an equation similar to equation (A. 10), new [E] 
and {F> transfer matrices like those in equations (A. 11) and (A. 12) 
can be defined as follows 

[Ej] = [I] +£l ([A,*] + 4[A] + [A i+1 ]) 


+ 



([A] CAiD + [A] [A] + [A 1+1 ] [A]) 


+ ^ (DO [A] [A-j] + [A i+1 ] CA] [A]) 

hi 4 

+ ~& CA i+l ] m CA i ] (4.12) 


{FiJ = 111 ({B-f > + 4 {¥} + {Bj +1 } ) 

6 

h .2 

+ ([A] {B-j} + [A] m + [A 1+1 3 {F}) 

h,-3 _ _ 

+ _I_ ([A] [A] (Bj> + [A 1+1 ] [A] {B}) 

hi 4 

+ — ^ CA-j D [A] [A] (Bj } (4.13) 
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From these new definitions of the transfer matrices (4.12) and 
(4.13), the rest of the solution, equations (A. 13-17), can be directly 
applied as usual. 

In addition to solving linear problems, the Hunter method, using 
either a second-order or a fourth-order Runge-Kutta integration, can 
be used to solve nonlinear differential equations. This is 
accomplished simply by collecting the nonlinear terms in the [A] 
matrix (or the { B> vector) and iterating the solution until successive 
{ Y} values differ by no more than some predetermined (convergence 
criteria) small value. An initial guess can be provided or simply set 
to zero (solving the uncoupled problem), and each new set of { Y} 
values can then be used to modify the nonlinear [A] matrix (or the 
{B} vector). This method has been used to solve several nonlinear 
differential equations with a great deal of success, and achieving 
convergence in usually five or six iterations. 


4.2 Fourth-order Transfer-matrix Error Estimate 
Since the above solution method involves a fourth-order 
Runge-Kutta integration which is an application of a fourth-order 
Taylor series, the error induced in using the approximation can be 
estimated by evaluating the Taylor series remainder term [10] 


h-j 4 d 4 {Yj > 
4! dx 4 


(4.14) 


The fourth derivative of { Y) with respect to x can be evaluated by 
using the chain rule of calculus on equation (4.1). 
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CHAPTER 5 


COMPUTER PROGRAMS 

Two computer programs, written in the FORTRAN programing code on 
a personal computer, were developed to solve the equations developed 
in chapters 2 and 4. The first program, PROP, computes the cross- 
sectional properties for a thin open cross section; and the second 
program, SOLVE, solves the governing equations that are summarized in 
figure 3.1. Both of the programs were developed for the desk -top 
personal computer and were designed to run independent of each other; 
however, they are easily coupled through a data-base or spread-sheet 
program such as LOTUS-1-2-3. 

5.1 Program PROP 

Program PROP computes all of the required cross-sectional area 
properties that are described in chapter 2.4.2 by using the equations 
developed in appendix B. This program is written in FORTRAN-77 for 
the personal computer and can easily be up-loaded to a mainframe 
computer that uses a similar version of FORTRAN. The program only 
requires ordered pairs of coordinates that adequately define the upper 
and lower surfaces. These points need not be evenly spaced or have 
coincidence abscissa coordinates. 
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Figure 5.1 is a flowchart of the program PROP. The program, upon 
being started, reads from a user-defined file the upper and lower 
surface definitions. Next, the user is prompted to input the number of 
equally spaced integration points. This option is provided so that 
the user can determine if the numerical integration of the equations 
in appendix B have converged. The equal spacing of the integration 
points are determined by linearly interpolating between the input 
surface definitions. Thus, a reasonable description of the surfaces 
needs to be input. For the geometries considered in this analysis, 
the cross-sectional properties indicated convergences with 150 to 200 
integration points. After the number of integration points has been 
determined, the program then computes the chord length and the n 
equally spaced stations. From this data, the area and shear center 
are computed, then the cross-section properties are computed about the 
shear center. The user is then prompted to input an angle about which 
the properties are to be transformed. The transformed properties are 
computed and printed, and the user is prompted to input more data. If 
more data are to be analyzed, the program is reinitialized and 
restarted. 


5.2 Program Solve 

The Hunter transfer-matrix method described in appendix A and the 
fourth-order matrices developed in chapter 4 have also been programed 
in FORTRAN-77 for the personal computer. This program yields a solu- 
tion to the equation summarized in figure 3.1 using the boundary 
conditions in equations (3.17) and (3.18). 
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Figure 5.2 is a flowchart of the program SOLVE. This program 
requires as input the number of station and the properties shown in 
the matrices in figure 3.1 for each station that is defined. Upon 
starting the program, all inputs are read in a table format, and all 
undefined variables are set to zero. The program is initialized by 
setting the iteration counter to one and assuming all nonlinear 
variables are zero. This usually solves the linear uncoupled problem 
and is a reasonable estimate for the second iteration. The station 
matrices, [A] and (B), are assembled in a user written subroutine that 
has been precompiled with the main program. Upon assembling the [A] 
and { B} matrices, the other system matrices are assembled. Then using 
the boundary condition, the first station is related to the last, 
resulting in a single linear matrix equation that requires solving. 
Once the complete solution at station one is known, then equation 
(A. 11) is repeatedly applied until all of the (Y-j > ’ s have been found. 
Since the check values of the solution variables are preset to zero, 
the first convergence check always fails and the second iteration is 
automatically started. After the second iteration, subsequent itera- 
tions produce smaller and smaller errors until the convergence checks 
are satisfied. When convergence is achieved the solution is printed 
and the program terminated. 
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CHAPTER 6 


RESULTS OF APPLICATION 

To determine the effect of the centrifugal force stiffening on a 
rotating beam, two example cases are evaluated, each for a blade 
rotating at a constant angular speed of 475 revolutions per minute 
(RPM). The first case compares the transfer-matrix solution of the 
governing equations to a nonlinear finite-element analysis of a 
rotating 7- by 10-Foot Wind Tunnel fan blade that neglects the effects 
of the mass axis and elastic axis being eccentric. The second case 
includes the effect of the mass axis being eccentric from the elastic 
axis for the transfer-matrix solution and evaluates the influence of 
the nonlinear terms on the response. 

6.1 Nonuniform, Noneccentric, Twisted Fan Blade 

A special case of the fan blade geometry presented in appendix B 
where the mass axis eccentricity is neglected is analyzed using both 
finite elements and the transfer-matrix solution. Both models were 
analyzed using a set of steady aerodynamic design load that were 
developed by the Ames Research Center specifically for the 7- by 
10-Foot Wind Tunnel fan blades. 

6.1.1 Transfer-Matrix Solution 

The equations that are summarized in figure 3.1 were solved using 
the program that is discussed in section 5.2 and the cross-sectional 
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properties summarized in appendix B. The steady aerodynamic loads 
shown in table 6.1 were assumed to act the quarter-chord point; and, 
therefore, had to be transferred to the elastic axis which resulted in 
a distributed torque along the blade in addition to the transverse 
loadings. A converged solution, following the flow-chart shown on 
figure 5.2, was achieved in seven iterations. Since the steady 
aerodynamic loadings were discontinuous at station 96, a double 
station was used at this station. This was accomplished easily by 
utilizing the feature of reducing the transfer matrix to an identity 
matrix when hj = 0 ( or a double station). Taking advantage of this 
precludes the solution algorithm from tapering the applied loading to 
zero from station-to station. 


TABLE 6.1 - AERODYNAMIC DESIGN LOAD INTENSITIES 


STATION 

(in) 

AXIAL 

(lbs/in) 

TANGENTIAL 

(lbs/in) 

96 

34.55 

14.33 

106 

45.23 

15.94 

116 

55.08 

16.54 

126 

65.37 

17.09 

136 

74.78 

17.24 

146 

83.10 

16.93 

156 

88.97 

16.04 

166 

90.48 

14.44 

173 

89.91 

13.07 
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6.1.2 Nonlinear Finite-Element Solution 


As a means of performing an accuracy check on the transfer-matrix 
solution, a finite-element model of the fan blade using two-noded beam 
elements and the EAL [12] finite-element computer code was analyzed on 
the CYBER CY17-855 mainframe computer. The finite-element analysis 
assumes the reference axis of the beam is the centroid axis while the 
transfer-matrix solution uses the elastic axis as the reference axis. 
To avoid over complication for a test case, the mass axis off-sets 
were neglected. This allowed the same loads to be input for both 
analysis without additional transformations, with the exception of the 
torque. The transverse loadings were input to the finite-element code 
as distributed load; however, no such feature exists for moments. 

This difficulty was overcome at the expense of inducing a small error 
by inputing the distributed torque as a discrete torque at each grid 
point. A similar approximation was induced by discretizing the 
geometry into elements with constant properties (the transfer-matrix 
solution assumes variable properties between stations). After 
defining the model and loads, the geometric nonlinear analysis (GNA) 
runstream element of EAL was modified to allow for the previous itera- 
tion's displacements to be used in computing the centrifugal forces 
due to the fan blade rotating at 475 rpm (49.742 rad/sec). 

Introducing this additional nonlinearity, resulted in the EAL code 
requiring 31 iterations before convergence was achieved. 

6.1.3 Transfer-Matrix and Finite-Element Solution Comparison 

The solution results for the finite-element and transfer-matrix 
methods are compared. This method of validating the transfer-matrix 
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method for the rotating beam problem was selected due to the absence 
of any experimental or analytical data in the literature. All of the 
studies that were surveyed were centered around predicting and corre- 
lating the natural frequency response (frequency and mode shapes) of a 
rotating beam. If any displacements and loads were computed, they 
were for a nonrotating propeller or fan blade under static loads. 

The two solution methods could not be compared directly on all 
responses due to the basic formulation differences. Comparisons are 
made on the displacements, rotations, twist, and support reactions. 
Internal loads could not be compared due to the finite-element 
approach of discretizing variable properties. However, all of the 
other quantities were compared. Support reactions for both solution 
methods are shown in table 6.2; these results indicate that there is a 
excellent agreement between the two methods. The largest percentage 
difference was in the twisting moment (only 2.1 percent). This is 
expected since the pretwist of the blade is explicitly accounted for 
in the transfer-matrix method, and is only approximately accounted 
for in the finite-element solution by transforming the principal axes 
of the cross section for an element. The displacement results 
(figures 6.1 through 6.6) with the exception of the axial displace- 
ment, figure 6.1, of the blade, agree extremely well. What is seen 
here is the primary difficulty of using finite-elements to solve this 
type of problem. The application of the centrifugal force, using 
finite elements, results in the inertial force (a function of mass, 
location, and deformations) becoming a lumped nodal quantity which at 
the stress-free end of the blade (tip) is not stress free. Because of 
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this condition, the displacement quantities at the end of the blade 
using finite elements may be in error. 

TABLE 6.2 - SUPPORT REACTION COMPARISON 


REACTION 

TYPE 

FINITE ELEMENT 
RESULTS 

TRANSFER MATRIX 
RESULTS 

PERCENT 

DIFFERENCE 

T (lbs) 

173,205 

173,835 

- 0.4 

V y (lbs) 

- 1,455 

- 1,447 

0.5 

V z (lbs) 

5,382 

5,382 

- 0.0 

M x (in-lbs) 

33,911 

34,632 

- 2.1 

My (in-lbs) 

-367,055 

-372,021 

- 1.3 

M z (in-lbs) 

- 82,547 

- 82,980 

- 0.5 


With the noted exceptions, the finite-element and transfer-matrix 
solution method agreement is good; thus, this comparison is considered 
to verify the transfer-matrix solution of the governing differential 
equations for a rotating beam. 

6.2 Nonuniform, Eccentric, Twisted Fan Blade 
After verifying the solution method, the fan blade model that 
was used in section 6.1 was modified to include the offset of the mass 
axis from the elastic axis. The solutions for the fully coupled case 
(the case where the full set of equations as indicated in figure 3.1 
are utilized) and the partially coupled case (the case where all dis- 
placement nonlinearities are neglected) are compared. This comparison 
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illustrates the influence of the displacement nonlinearities that are 
indicated in figure 3.1. 

6.2.1 Fully Coupled Transfer-Matrix Solution 
The next solution that is performed is for the fully coupled, 
nonlinear, nonhomogeneous differential equations. This represents the 
design operating condition for the 7- by 10-Foot Wind Tunnel fan 
blades. The steady-state aerodynamic loads at the quarter-chord point 
are shown in table 6.1; these loads as indicated in section 6.1.1 are 
transformed to the elastic axis. By including all of the displacement 
nonlinearities and the mass axis eccentricity, the solution converged 
in 12 iterations. 

6.2.2 Partially Coupled Transfer-Matrix Solution 
To determine the influence of the nonlinear terms shown in 
figure 3.1, all of the nonlinear terms that involve displacement and 
rotation terms are neglected; this leaves only the centrifugal axial 
force, T, coupling term in the bending moment equations. Neglecting 
the nonlinear displacement terms, the solution converged in three 
iterations. This is easy to see by noting that the axial force, T, 
equation is weakly coupled to the other equation through the axial 
displacement. The axial displacement term is usually neglected [2], 
and, therefore should have little influence on the net axial force 
solution. Thus, the axial force solution is nearly exact as a result 
of the first iteration; and upon using the first iteration's results 
in the second (neglecting the nonlinear displacement terms) the second 
iteration should result in a reasonable accurate solution. The third 
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iteration is used to insure that the absolute value of the difference 
between successive iterations is less than the required criteria. 

6.2.3 Fully Coupled and Partially Coupled Transfer-Matrix 

Solution Comparisons 

Since the number of interations that were required to achieve a 
converged solution were found to vary drastically when the displace- 
ment nonlinear-terms were considered, a comparison of the solutions is 
made. The more notable effects of these nonlinear terms were observed 
to influence the displacements and rotations, and had little or no 
effect on the forces and moments. 

A comparison of the axial displacement is shown in figure 6.7. 

The inclusion of the nonlinear displacement terms has a stiffening 
effect that reduce the tip displacement by approximately 25 percent; 
this is due to the blade not untwisting as much when the nonlinear 
terms are included. 

By considering the nonlinear displacement effects on the 
transverse deflections (figures 6.8 and 6.9), the deflections 
increased on the order of 3 percent. This indicated an insignificant 
softening of the fan blade. 

The most significant effect on the response is noted when the 
cross-sectional twisting of the fan blade is evaluated. This effect 
is essentially an order of magnitude stiffer, see figure 6.10. It is 
probably for this reason alone that these terms are necessary to 
include in performing a loads analysis of a rotating fan blade. For 
this analysis, the aerodynamic loadings were considered to remain 
constant and independent of the deformed shape of the blade. 
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Normally, however, this is not the case, since the aerodynamic 
loadings are generally a function of the angle-of-attack of the blade 
(relative to the free-stream) . The lift and drag loads are then 
displacement dependent. If the aerodynamic loads are a function of 
the blade's displacements, then this condition can be included by 
introducing an additional nonlinear term in the load vector which is 
iterated along with the other terms or by reformulating the problem as 
a beam on an elastic foundation and including the displacement 
dependent loads as linear terms in the A-matrix. 

Like the transverse displacements, the bending rotations, by 
including the nonlinear terms, demonstrates a softer response, 
see figures 6.11 and 6.12 for this effect. These terms are strongly 
coupled in the moment equations to the mass axis eccentricities, and 
therefore would be expected to influence the solution. 

Unlike the displacement results, the internal equilibrium loads 
are not strongly affected by the presence of the nonlinear displace- 
ment term. Figures 6.13 through 6.18 show little or no change in the 
load results. In fact with the exception of the transverse shear, Vy, 
and beam bending moment, M z , (figures 6.14 and 6.18) plots of the two 
solutions completely overlay each other. These two loads are coupled 
and influenced by the y-inertial force component which, according to 
equation 2.43, is dependent upon the v transverse displacement and 
accounts for this small variation. 
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CHAPTER 7 


CONCLUDING REMARKS 

Nonlinear equations of motion for the coupled elastic bending and 
torsion of twisted nonuniform rotating beams are derived. In addi- 
tion, a transfer-matrix solution method to solve these nonlinear 
equations is developed and presented. Two case were evaluated. The 
first case compares the transfer-matrix solution method to results 
obtained from a nonlinear finite-element code. 

The nonlinear equations were developed by neglecting all but the 
first-order terms for the strain, and retaining all other nonlinear 
terms in the derivation. The resulting equations are for the coupled 
bending-torsion steady-state response of beams rotating at a constant 
angular velocity in a fixed plane. 

The modified Hunter transfer-matrix method was verified by 
comparing results with solutions from a geometric nonlinear finite- 
element computer code. These results were shown to agree quite well, 
and that as is the case with any special purpose analysis when com- 
pared with a general purpose code, the transfer matrix yielded a more 
accurate representation and solution of the problem. An analysis of a 
proposed new fan blade design for the 7- by 10-Foot Wind Tunnel at the 
Langley Research Center was performed, and the effect of including 
nonlinear displacement terms on this solution addressed. 
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As a result of performing this analysis, two computer codes were 
developed for a desk top personal computer. The first program is a 
general purpose two-point boundary value problem solver that was used 
to solve the nonlinear, coupled governing differential equations of 
motion for a rotating beam. The second program performs a numerical 
integration for the cross-section properties for a thin open cross 
section that can be defined by an upper and lower surface. 

This thesis develops the effective and efficient analysis 
techniques that consider the unique features common to rotating 
systems to determine the steady-state response of rotating blades. 
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APPENDIX A 


THEORETICAL DERIVATION OF THE HUNTER TRANSFER MATRIX METHOD 

The method that is presented here is a solution procedure that 
was developed by Dr. William F. Hunter, and used by him for several 
years to solve various homogeneous and nonhomogeneous two-point 
boundary-value problems. Because to date the method has not been 
documented, it is developed in detail in this appendix and referred to 
as the Hunter method. 

The development of the transfer matrix for the linear coupled 
matrix equation of the form 

= [A] { Y> + {B} (A. 1) 

dx 

is described. By expressing the equations in matrix notation, 
utilizing the transfer matrix as an operator, and applying the 
boundary conditions, the coupled linear differential equations are 
solved completely in two passes over the solution domain. The first 
pass is to satisfy all of the boundary condition (natural and 
geometric) at the first station; and the second pass develops the 
solution at each station within the problem domain. 
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Consider a system of linear nonhomogeneous first-order 
differential equations that can be expressed in matrix notation as 


d - { - Y i x) - ? = [A(x)] { Y( x ) } + {B(x ) } UK. 2) 

dx 


defined on the closed interval x^<x<x n . 

Using subscripts to denote stations and primes to denote 
differentiation with respect to x, equation (A. 2) at station 
i (x = x-}) is written as 


{Yi 1 } = 

CAi] 

{Yi} + {Bi } 

(A. 3) 

The boundary conditions at the 

end-| 

points (x = xi and x = x n ) 

are al 

expressed in matrix notation as 



[C] 

{Yl> 

= {P} 

(A. 4) 

CD] 


« {Q} 

(A. 5) 


Approximating the value of { Y> at the station (i+1) as follows 


where 


{Y i+ 1> = { Y-j > +ihi ({Y i '} + {Y 1+1 ‘}) 
2 


hi ~ x i+l " x i 


(A. 6) 


(A. 7) 


then equation (A. 2) may be written at station (i+1) and substituted 
along with equation (A. 3), into equation (A. 6) to give 


Now, the Y-j+j term on the right-hand side of equation (A. 8) may be 
approximated by using a Taylor series expansion about Yj . Retaining 
only those terms up through the first derivative, Yj+j becomes 

{Y i+ i} = (Y-j } + {Yj 1 } 


or 

{ Y i + 1 } = {Y-j} + ([A-j] {Y-j } + (B i }> (A. 9) 

Substituting equation (A. 9) into equation (A. 8) gives 

h ■ u .2 

(Y i+ i) = {[I] + j- ([A-j] + [A-j +1 ] ) + -I- [A-j +1 ] [Aj]) {Y-j} 
h* h 

+ Y ( {B-j } + {B i+1 }) + ^- [A i+1 ] {B-j } (A. 10) 

where [I] is the identity matrix. It can be shown that the above 
result is a second-order Runge-Kutta integration (see ref. 8) for a 
system of equations expressed in matrix form. 

Equation (A. 10) may be rewritten as 

{Y i+ i> = [Ei 1 {Y-j } + {Fj } (A. 11) 

where 

2 

[ Ei ] = [I] + ([A-j] + [A-j+il ) + [A-j +1 ] [A-j] (A. 12) 

h * 

{Fj} ({B 1 } + {B i+1 }) + -y- [Aj +1 ] {Bj } (A. 13) 


50 



The matrix [E-f] is known as a transfer matrix, since it relates 
the values of the state variables at station (i+1) to those at the 
station i The transfer-matrix approach is often used to determine 
the natural vibration characteristics of beams However, for a 
nonhomogeneous system of equations, the solution is complicated by the 
appearance of the matrix {F-j} in equation (A 11) Because of {F-j}, an 
approach much different from that of the natural vibration (homogen- 
eous) problem is needed to relate the conditions at one boundary to 
those at the other boundary It is at this point that the shooting 
methods will iterate the equations until conditions at both boundaries 
are satisfied The Hunter transfer-matrix approach as outlined here 
eliminates this difficulty by employing a systematic application of 
equation (A 11) from one boundary to the other; thus, the boundary 
conditions are satisfied directly The development is easily seen by 
using equation (A 11) and expanding it out for a few intermediate 
stations between boundaries For example, 

{Y 2 > = [Ej] {Y 2 > + {Fj} 

{Y 3 ) = [E 2 ] {Y 2 } + {F 2 } = [E 2 ] ( [E i 3 { Y x } + {F X >) + {F 2 } 

{Y 4 > = [E 3 ] { Y 3 } + {F 3 } = [E 3 ] ([E 2 ] ( [Ej ] { Y 2 > + {Fj}) + {F 2 }) + {F 3 } 

Now, the above can be refactored, and {Y^} can be expressed in general 

as a function of {Yj} For example, the expression for {Y 4 } above 
becomes 

{Y 4 } = [E 3 ] [E 2 ] [Ej] {Yj} + [E 3 ] ([E 2 ] {Fj} + (F 2 } ) + {F 3 } 
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In general, the relationship can be expressed as 


{Y*} = [Gj] {Y]l> + {H i } (1-1,2 n) (A. 14) 

where 

[G-j ] = CE^.i] [G 1 (i=2,3,...,n) 

CGiD « Cl] (A- 15) 

{H n *> = [E^] {H 1 _ 1 } + £F 1 _ 1 > (i=2,3, . . . ,n) 

{Hi> ■ 0 (A. 16) 

With equation (A. 14), the state vector at one boundary can be related 
to the vector at the other boundary. Thus, equation (A. 15) for i=n 
becomes 


{ Y n } = CG n ] { Yj } + {H n } (A. 17) 

The boundary conditions must be applied in order to determine { Yj > - 
Equation (A. 17) may be substituted into equation (A. 5), and the 
results combined with equation (A. 4) to give 


~ CC] - 

{ Yi } = 

{P} 

ID] [G n ]_ 

IQ} - [D] { H n > _ 



(A. 18) 


Since the coefficient matrix is nonsingular, equation (A. 18) can be 
solved for { Yi> - The solution for { Y-j } ( i >1 ) is obtained simply by 
applying equation (A. 11) repeatedly. 

It is worth emphasizing that the Hunter method, given above, does 
not require any iteration, and obtains the solution in a direct 
manner. As with other transfer matrix methods, very little computer 
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storage is needed. The method requires stepping through the problem 
domain from one boundary to the other twice. The first step-through 
is made to obtain (Yi> , and the second step-through yields the 
solution. 

Another feature is that the integration process conveniently 
handles any discontinuities in the physical properties (such as beam 
mass or stiffness, for example) of the problem by allowing double 
stations at any x-,- (note that the transfer matrix [E-j] reduces to the 
identity matrix for h-j = 0, see equation (A. 12).) 
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APPENDIX B 


CALCULATION OF 7- BY lO-FOOT WIND TUNNEL 
FAN BLADE SECTION PROPERTIES 

In order to perforin the analysis of the 7- by 10-Foot Wind Tunnel 
fan blade new design, the cross-sectional properties as defined by 
equations (2.65-74) need to be calculated about the shear center of a 
cross section relative to the beam's y-z coordinate system. (See 
figure 2.1.) Since the definition of a fan blade is usually given by 
coordinates that define the upper and lower surfaces, equations 
(2.65-74) can be rewritten in a form that will allow for numerical 
integration of the equations. The purpose of this appendix is to 
present the modified and transformed cross-section integral equations 
as well as the numerical results when applied to a new 7- by 10-Foot 
Wind Tunnel fan blade design. 

A typical cross section of an ai rfoil -shaped fan blade and 
associated coordinate systems are shown in figure B.l. The cross- 
sectional coordinates of the lower ( x ^ ) and upper (x u ) surfaces are 
given in the ip - x system. When the shear center is located, the 
yi - z\ system is used to compute the section properties parallel to 
the ip — X system but about the shear center (ip s , \ s ). Since a typi- 
cal cross section can be twisted about the shear center relative to 
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the y-z system, the cross-sectional properties need to be transformed. 
The section properties transformed to the y-z system (shown in 
figures 2.1 and B.l) are then calculated. 

Blade description data is assumed to be defined in the ip — X 
system. The area and shear center location are computed in this 
system by the following equations. 


c 

A = J* (A u - \ t ) d * 
0 


(B.l) 


♦s 


£ ip( x u - d 

o 

c 

S (Ay - A^)^ (j \|) 


(B.2) 


_ (A u - X t ) 


X S ~ 2 “ 


(B.3) 


4’ = i's 


where c is the chord length. 


Equation (B.2) is developed by considering equilibrium of a 
differential element under bending and requiring that the sum of the 
internal moments about the point called the shear center to vanish. 
Equation B.3 is an approximation that is used in this analysis for a 
thin airfoil shape. 
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The other section integrals are 


y i =7 j U - ^ ^ x u - x ^.) d ^ 
A o 


(B.4) 


- _ i r /u 2 x x 

Z1 - 7 J (— ’ Xs l 
C 0 


Y + x s x * } d * 


(B. 5) 


Iy 1 yi -\ S C(X U - x s )3 - (X* - x s) 3 3 d * 
o 

c 

I Z1 Zl = / (* - ts * 2 < x u - x *> d ♦ 

0 


(\|> - i|> s ) CU u - X s )2 - (x A - X s ) 2 ] d ^ 

(B.9) 

( B - 10 ) 


l yi zi 


■s 


Pi = yi a 
•P 2 = Z! A 


(B.6) 


(B.7) 
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^3 -/ {C(X U - X s > 2 - (X t - X $ )2] (l^) 2 


+ - C(X U - X S )4 - (\ z - X s )4] } d rp (B.ll) 


c 

p 4 ~S ~ ts^ 3 ( x u " + ( — 2 ^ x u " x s^ 3 

o 

- (X A - X s )3]} d \p 


c 

^5 ~S {(♦ - ( x u - X fc) + [( x u • x s> 5 - ( x £ 
0 


(B. 12) 


X s )5] 


— j [(X u - X s )3 . (x £ - X s )3]} d i> 


(B. 13) 


The torsional stiffness constant [11] can be calculated by 


J = 


F 

1 


1 + 


4 F 

3 A c 2 


where 


F =/ ( x u ~ X £ ) 3 dc (B.14) 


The following transformation rule relates the y\ - z\ system to 
the y-z system. 
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y = yi cos (e t ) + zj sin (B t ) 

z = -yi sin (6-^) + z\ cos (B^ (B.15) 

Using the transformation indicated in equation (B.15), the 
cross-section properties (equations (B .4-13) ) in the y-z system 
become 

y m = 7i cos (B t ) + zi cos (8 t ) (B.16) 

z m = - 7i sin (B t ) +1! cos (3 t ) (B.17) 

lyy = I yl yj COS 2 (g t ) + Izi Zj Sin 2 ( 0 t ) 

- 2 I yl Zl cos (B t ) sin (B t ) (B.18) 

I Z z = r yi yi sin 2 W + l z\ z x c °s 2 (& t ) 

+ 2 I yi Zl cos (Bt) sin ( B-t ) (B.19) 

I yz = I yi yi [cos (Bt) sin (B t )] - I Z1 Zl [cos (&t) sin (B t )l 

+ I yl Zl [cos 2 (Bt) " sin 2 (e t )3 (B.20) 

Pi = y m A (B.21) 

P 2 = z m A (B.22) 

p 3 = - ^4 sin (B t ) + P 3 cos (B t ) (B.23) 

P 4 =^4 cos (B t ) + "P 3 sin (B t ) (B.24) 
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The properties. A, J, and P 5 , are invarient, and, as such, they are 
not effected by coordinate transformations. 

A sketch of a 7- by 10-Foot Wind Tunnel blade is shown in figure 
B.2. The stations indicated along the span of the blade are the 
points where the cross sections were cut and properties calculated. 
Figures B . 3-18 show the cross section and indicates the cross- 
sectional properties in the y-z systems. For this analysis the fan 
blade from station 57.375 down toward the center of rotation was 
considered fixed due to the rigid attachment to two thick steel disks. 
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(a) Bending in the x-y plane 



(c) Extension of elastic axis 


Fig 2.3 Cross-sectional deformations 
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(a) x-y plane 



Elastic axis 



y - z 0 
J m m 

i + y m 0 
m J m 


Cross-section mass center 



(P + O) 


(c) y-z plane 


Fig 2.6 Equilibrium of differential bending element 
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65 


Fig 2.8 Contribution of longitudinal stress 
to the internal elastic torque 






= [A]| Y I ♦ IBI 
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Fig 3.1 Governing differential equations in matrix format 



('Read upper/lower surface coordinates 


c 


E 


Input number of integration points 
n 

I 


Compute chord length and n equally spaced stations 

1 — — 


Compute area and shear center 

i 


Compute cross-section properties about shear center 

l 


Input transformation angle 
l — 


Recompute cross-section properties 
about transformed system 


i 


/ Output computed properties/ 


More 

"section property calculations required 
? 


Yes 


No 


Fig 5.1 Flow chart of program PROP 
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Fig 5.2 Flow chart of program SOLVE 







Axial, u, displacement 



Beamwise, v, displacement 



Fig 6.2 Comparison of beamwise displacement 
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Chordwise, w, displacement 



Twist, <t>, rotation 



Fig 6.4 Comparison of cross-sectional twist 
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Bending, 0 y , rotation 



O Transfer matrix 
□ Finite elements 


Fig 6.5 Comparison of rotations 


Rotations, 


.05 
.04 h 
.03 )- 
.02 
.01 
0 
01 
02 

03 

04 


Bending, 0 Z , rotation 


O Transfer matrix 
□ Finite elements 



Radial, u, displacement 



Beamwise bending, v, displacement 
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Chordwise bending, w, displacement 



Twist, O, rotation 
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Rotation, 

radians 


Rotation, 

radians 


Bending, 6 , rotation 



Bending, 6 , rotation 
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Axial, T, force 


O w/o nonlinear terms 
□ w/ nonlinear terms 



Beam station, in. 
Fig 6.13 Axial force, T 


Shear, V , force 
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Moment, 

in-kips 


Moment, 

in-kips 
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Fig B.3 7- by 10-Foot Wind Tunnel fan blade station 57.375 
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Fig B.4 7- by 10-Foot Wind Tunnel fan blade station 60.000 
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